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. . . The Kronecker channel model of wireless communication is analyzed using statistical mechanics 

0^ ' methods. In the model, spatial proximities among transmission/reception antennas are taken into 

^•ZJ , account as certain correlation matrices, which generally yield non-trivial dependence among symbols 

^<ZJ • to be estimated. This prevents accurate assessment of the communication performance by naively 

using a previously developed analytical scheme based on a matrix integration formula. In order to 
resolve this difficulty, we develop a formalism that can formally handle the correlations in Kronecker 
^ ■ models based on the known scheme. Unfortunately, direct application of the developed scheme is, 

' in general, practically difficult. However, the formalism is still useful, indicating that the effect 

of the correlations generally increase after the fourth order with respect to correlation strength. 
Therefore, the known analytical scheme offers a good approximation in performance evaluation 
when the correlation strength is sufficiently small. For a class of specific correlation, we show that 
the performance analysis can be mapped to the problem of one-dimensional spin systems in random 
fields, which can be investigated without approximation by the belief propagation algorithm. 
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I. INTRODUCTION 



>. 

04 ^ Recently, in the field of information science, techniques for efficiently handling systems with large amounts of data 
OO . are strongly required, and statistical mechanics have attracted a great deal of attention. The number of applications 
in information science to which analytical schemes in statistical mechanics can be applied is increasing, and such 
, applications offer a variety of consequences [l[ , some of which are not possible by standard techniques in information 
0^ ' science. Information processing is a notable example. 

, In the present study, we investigate wireless communication systems. Multiple - input multiple - output (MIMO) 
systems and code division multiple access (CDMA) systems in wireless communication have mathematical structures 
that are similar to those of disordered spin systems in physics, and analytical tools in statistical mechanics, such 
as the replica method and mean field approximations, have enabled performance analysis and improved processing 
algorithms for actual communication systems 0, [li 0i [1, fl 0, these studies, the communication process is 

described by a linear equation with transmitted signals h G and received signals r G using an. L x K channel 



matrix H = (Hik) G C'^'^'^ and noise rj E as 

r = Hb + arj, (1) 

where describes the noise power. Throughout the present paper, matrices and vectors are denoted in bold. 
In the above equation, the dimension K represents the number of multiple transmission antennas in the MIMO 
system, whereas L corresponds to the number of reception antennas. Clarifying the feature of the above-mentioned 
communication channel by the standard method of information theory is technically difficult because of the randomness 
in H and the discreteness of variable b. However, statistical mechanical analysis enables us to avoid such difficulties 
in the limit of infinite system size. 

In previous studies based on statistical mechanical methods, channel matrix H was characterized by a property 
that the cross correlation H^H can be handled as a typical sample from a rotationally invariant matrix ensemble, as 
follows: 

H'^H = UDU\ (2) 



where C/ is a sample randomly chosen from the uniform distribution of iiT-dimensional unitary matrices, and £> is a 
if-dimensional diagonal matrix. If D has an asymptotic and deterministic eigenvalue distribution pr>(A) with a large 



2 



matrix size limit L,K ^ oo while keeping f3 = K/L finite, the features of this channel can be characterized by p(A), 
in conjunction with the replica method, and the performance of the channel can be assessed [Tol . [Tl| by a matrix 
integration formula [13, [13 H, [IB] , which is defined for p(A). 

However, one problem remains. For the simplest case in which each element of random matrix H is drawn from an 
independent and identically-distributed (i.i.d.) Gaussian distribution, the property of rotational invariance concerning 
the cross correlation is satisfied. However, such property does not necessarily hold for general matrix ensembles of 
MIMO systems. For instance, in the Kronecker model [l^, which is one of the standard models in the theory of 
wireless communication, the elements of the channel matrix are not drawn from an i.i.d Gaussian distribution, but 
are instead drawn from an L x _ftr-dimensional joint Gaussian distribution. More precisely, the channel matrix H is 
described as 

H = y^R^Sy^Rt, (3) 

where each component of an L x ii' rectangular matrix H — (Sj^) is drawn from a complex i.i.d. Gaussian distribution: 
P(Eik) = L7r-ie-^l^"=l^(l < I < L,l < k < K). € C^^ and Rt € C^' are L- and i^-dimensional determin- 
istic matrices, which are Hermitian and indicate correlations among reception antennas and transmission antennas, 
respectively. In a previous studv[loj. we analyzed this system by means of the matrix integration formula. However, 
for this system, the matrix ensemble is not rotationally invariant and, accordingly, the result of performance analysis 
via the matrix integration formula may not hold exactly. 

One of the goals of the present paper is to develop a scheme that can handle the dependence on y'R^ and \/Rt in Eq. 
explicitly. In other words, the method developed herein relies on the direct integration of each matrix element in 
Gaussian random matrix H. The results of analysis for mutual information indicate that the roles of the deterministic 
matrices \/Rr and \/Rt are different. As will be shown later herein, the dependence on ^/R^ can be treated using 
the matrix integration technique, whereas the dependence on ^/Rt must be handled more carefully. The developed 
scheme can also be used to construct a practical demodulation algorithm. Another goal is to compare the performance 
of the Kronecker channel Q via a novel analysis with the performance of the matrix integration formula applied to 
the entire cross-correlation matrix H^H, as we demonstrated in Reference [lo| . The two formulations are found 
to yield different result, which means that the application of the matrix integration formula to the cross-correlation 
matrix H , in general, does not yield correct results. However, when correlation among transmission antennas, or 
when the off-diagonal element of the deterministic matrix ^/Ri is sufficiently small, discrepancy between results of 
the scheme developed herein and that based on the matrix integration formula increases only after the fourth order 
with respect to correlation strength, implying that the formulation based on entire matrix integration yields good 
approximate results. This suggests that although the matrix-integration technique is generally an approximation, 
this technique is practically useful when the correlation is small, because the matrix integration method enables the 
system to be characterized using only a few macroscopic variables, which significantly reduces the computational cost 
for analysis. 

The remainder of the present paper is organized as follows. In Section |TT1 we provide basic tools for performance 
analysis and propose a novel approach to analyze the Kronecker channel model. The analytical results differ from 
those obtained by matrix integration. In addition, we compare two results by a method of perturbative expansion 
with respect to the correlation parameter, and a discrepancy appears in the fourth-order coefficient of the correlation 
parameter, which indicates that the discrepancy is small when the correlation is small. In Section lllll we show that the 
demodulation algorithm can be constructed from the minimization scheme of Gibbs free energy without knowledge of 
the matrix integration. We present the experimental results of the demodulation algorithm for the Kronecker channel 
in Section HVl As a special case, we consider a system with a tridiagonal form of Rt, where the analytical scheme can 
be used for the random-field Ising chain. The results of a numerical experiment confirm the validity and usefulness 
of the proposed scheme. The final section presents a summary of the present paper. 

II. ANALYSIS 

Let us start with the communication channels described by Eq. ([T]). For the noise, we assume that rj is drawn 
from a white normal complex Gaussian distribution P(rj) = Tr^^e^'*'! . Each component of the iiT-dimensional 
transmit vector b is generated from an i.i.d. information source and modulated. For simplicity, the modulated 
components 5^ are quantized to one of the elements in a set B. For instance, for S'-phasc shift keying modulation 
B = {e^^^'^Z ^}{s = 0, 1, S — 1). As special and well-known cases, B = {±1} for binary-phase shift keying (BPSK) 
modulation, and B = {±l/\/2 ± i/\/2} for quadrature-phase shift keying (QPSK) modulation. The prior of the 
transmit vector is denoted by P{b) = Yl^=i Pibk)- Here P(6fe) = l/I^Bj, and \B\ is the number of elements in B. 
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As mentioned in the introduction, we investigate the Kronecker model described by the matrix of Eq. ([3]) . In order 
to apply statistical mechanical schemes to the analysis of communication systems, we allow the number of antennas L 
and K to be sufficiently large while keeping (3 — K/ L finite. Next, let us assume that for the matrices Hr and Rt that 
there exist deterministic distributions pn^X) and pr^{X), respectively, in the limit of infinite number of antennas. In 
addition, we assume that both distributions have compact supports and finite moments, which affects the applicability 
of the matrix integration formula. 

In the following, we consider only the case in which the receivers know the channel matrix H and the noise power cr^ 
in advance. The performance of the communication channels can be analyzed by estimating the mutual information 
between transmitted signals b and the received signals r, denoted by Th- For MIMO systems, we have 



Jh ^ f drZ(r) In Z(r)-iln(7ra2)-i, 

K JcL p p 

where Z{r) = TrP(6)- 



exp 



Hb\ 



(4) 



In this article we use nat unit for mutual information and entropy. In statistical mechanics, Z{r) serves as a partition 
function, which depends on quenched randomness if, and Xfj is considered to represent the free energy. 

Following the standard technique, we use the replica method to take the average over the channel matrix H in the 
mutual information: 



Ih ^ - lim -^^In 

n^o On K 



(5) 



where • • • denotes averaging over the distribution of channel matrix H . For n = 0, 1, 2, 



we have 



/ drZ"+i(r) = / drfr fTrP(6'^)^4\^exp 



(7rcr2) 



2\L 



exp 



Tr (H^'HL] 



(6) 



where the fC-dimensional square matrix L is defined by Lkk' = ELo KK* " Eafc K^t'/i^ + !)■ This can be rewritten 
as Lkk' = J^ab^k'P'^^^k* by introducing the (n + l)-dimensional projection matrix 7"^'' = 5'^'' — l/(n + 1). Substituting 
channel matrix H, given in Eq. ([3]), for the Kronecker model and integrating with respect to H, we obtain 

^^drZ«+i(r) oc |^J]TrP(6'^)^dct(^/iA' + ^fir®/R^iV^) 

= Pib")^ exp L J dXpR^ (A) Tr In (^Ik + ^ ^/RtLy/Rt 



Va=0 

J dQexp XTrGstfi.H (^^^q) +lnn(")(Q) 



RtL\/ Ri 



(7) 



where Id is the ZJ-dimensional identity matrix, Q is an (n+ f )-dimcnsional matrix, and ® represents the Kronecker or 
direct product. Note that the trace in Eq. ([7]) is if-dimensional in the second and third lines and is (n+ l)-dimensional 
in the last line. In the equation above, the following functions are defined: 



n(") (Q) = <^ n P^^') \ n Sib^^Rtb" - KQaa) \ { n Sib'^^Rtb' - KQab) 
GsfR^siA) = dXpRAX)\n{I-(3XA), 



(8) 



(9) 



where pr^{X) in the function G^fR^s is the eigenvalue distribution for the matrix R,-. The function GntH^s is the 
function obtained from the matrix integration formula over unitary matrix Haar measure 



J dUexp (TrHtiJrHA) ~ exp (if Tr Gstfl^s (A/K)) . 



(10) 
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Here, unitary matrix U to be integrated is defined as a JC-dimensional unitary matrix that diagonalizes the random 
matrix product as Ri-S — DU, where D is a diagonal matrix, and A is an arbitrary iiT-dimensional matrix. 
The result, given by Eq. ([7]), indicates that the entire set of unitary matrices that appear in the diagonalization 
of all possible random matrix products R^S coincides with the entire set of unitary matrices, which guarantees 
that the matrix integration formula over the unitary matrix is applicable only to the H^Hi-H part of the entire cross- 
correlation matrix H. This is because the multiplication of the matrix H serves as a unitary transformation. Note 
that we cannot apply the same argument to the entire cross-correlation matrix H = \/ RtS'^ RyS^/ Rt because the 
transmitter correlation matrix Rt breaks rotational invariance, and more careful treatment is required, as described 
below. 

By using the saddle point method and assuming replica symmetry as q = Qat for a ^ b and q + x — Qaai we can 
evaluate the replicated partition function in Eq. ([7]) after introducing auxiliary variables q + x ^-iid — 2g for the delta 
functions of the diagonal and off-diagonal matrix elements, respectively. 



- In n(") (Q) = Extr <^ nxx + {x~ {n + l)q){x+ {n + l)q) 

-K q.X 



n 
ah 



(11) 



The saddle point condition for q yields q = x/('^ + 1)- After performing Hubbard-Stratonovich transformation, 

K 



exp 



^-^ n + 1 ^-^ 



a=0 



{n -f l)x 



dr' exp 



^Rtb"" 



0=0 



(12) 



we have 



^lnn(")(Q) = Extr 



nxx 



K 



In / dr' I TrP(6) 



K 



exp 



'Rib 



n+l 



nln 



(13) 



Combining this equation with the remainder of the replicated partition function and noting that the matrix VQ has 
a single zero eigenvalue and n-degenerate x under the replica-symmetric condition, we obtain the final expression of 
the mutual information, as follows: 



= Extr(-G3tfl,3(-4) -T-^l^n("^(^) 
x-x \ \ / on K 

A 



x-x L 

= Extr(G3tH..HW + ^Rt 



n=0 



(14) 



where GsfR.s (A) is the Legendre transform of Gs^r^sWi GntHrS W = Extr^ {Ax — G'stH^H(x)}i and Ir^ (x) is 
defined as 



iRt (X) = -J^ 



dr' i Tr 

b 



Pib) (J) 



K 



exp 



'Rtb 



ln<^TrP(6)^ exp 



K 



'Rtb 



-ln( - 1 -1. 



(15) 



Equation (jl4|) is the primary result of the present paper. As we demonstrate in Section IIVI it is convenient to use 
iRt ix) foi' the discussion of the performance of the channel. 

In the following, we consider three items. First, Eq. (|14p provides a physical meaning for the performance analysis 
of the Kronecker channel. The term Ir^ of the right-hand side corresponds to the mutual information of a channel 
after rescaling of r'. 



(VA/a) 



(16) 



(See Eq. (U)) where rj' is a iiT-dimensional normal complex Gaussian noise. Here, we let A be a random variable 
that obeys probability distribution P(A) ~ exp KG^iR^sW 



Then, Eq. (|T4|) means that, in the limit, K 
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oo Th corresponds to the average of the exponential of the mutual information, exp(_ftr/jij (A/cr^)) over A. The 
extremization of Eq. p4|) implies that the balance of the two A-dependent functions G'stRjH(A) and /^^(A/tT^), which 
are dependent on correlations among reception antennas and among transmit antennas, respectively, is significant in 
the determination of Th ■ 

Second, Ib^ (x) can be evaluated using the following approximation method. After performing unitary transforma- 
tion of the matrix Rt to U'^ RtU, we take the average of lut Rtuix) over unitary matrix U (denoted by lu"* Rtuix) 
in the following) as in the case for the matrix Hi2iH^. In a manner similar to the evaluation of Th, we have 

IWR^uix) = Extr{GR,(A)+/,(Ax)}, (17) 

where Ii{x) is the mutual information of Eq. (jlSp after the substitution of Rt — I, which can be decomposed to the 
mutual informations of multiple one-dimensional channels. 

Third, if the correlation among transmission antennas is sufficiently small, we can perform a perturbative expansion 
of Ir^ . Let us consider the case in which the matrix Rt is expressed as Rt = I + pR with a real small parameter p 
and i^T-dimensional matrix ii, the diagonal elements of which arc all zero. After expansion, we have 

Ii+pr(X) = Ii[X) ^ J^UAX)} +— j^Uiix)} +•••, (18) 

where Ij{x) is the derivative of Ii{x) with respect to x- Similarly, expanding the approximate mutual information 
h+pU'' Ruix), we obtain the same result up to the third order of p. However, a discrepancy appears starting from 
the fourth-order coefficient. In the case of QPSK modulation, this discrepancy is expressed as (see the Appendix) 

Ii+pr{x)-Ii+pWru{x) = -^;^E{(^')*^-^^}V^/(x)-^;(x)'}«(x)}' 

||5:{Re(R.,)^ + Im(i?,,)4}j (^|_^;(^)_^;(^)2}2^£(|)!^ +0(^5)^ (19) 

where C(x) is a function that depends on P{b) as well as x- This indicates that the approximate evaluation of 
mutual information by matrix integration yields a good result if the perturbation parameter p is sufficiently small. 
Under this condition, the evaluation using matrix integration described in [lo| has an advantage in that it provides a 
good approximate solution that is more convenient than the exact evaluation of mutual information for the channel 
r = Hb + r). 

As described in the Appendix, we can prove that {— // (x) ~ -^/(x)^} ^ Oi and, accordingly, the right-hand side of 
Eq. pop becomes non-positive for a wide class of P{b), including QPSK modulation, which means that approximate 
evaluation by Ii+pu'' Ruix) gives an upper bound of Ii+pr{x) up to the forth order of the correlation parameter p. 

III. DEMODULATION ALGORITHM 

For practical communication, it is also significant to construct a computationally feasible demodulation algorithm. 
For inference of original signal b from received signal r and channel matrix H, it is necessary to evaluate the following 
quantity: 

m = E bP{b\r,H). (20) 

However, it is computationally difficult to numerically evaluate m from this expression. Key to the practical solution 
of this problem is the use of the Gibbs free energy for the communication channel: 



$(m) ^ Extr I - In Tr P(b|r , H) exp[Re(/i^ (b - m))] | 



(21) 



and the quantity m can be estimated as the argument of the extremized Gibbs free energy. Substituting P(b|r, H) = 
P{b) exp [— |r — ffbp/fj^] /Z with Z being the normalization and H = A/R^SV^Rt we have 



|r Hm 



2 



$(m) = Extr <j ^ — - In Tr ( P{b) exp 



^R^S\/R^{b — m 



|2 



Re{h\b-m)) 



InZ. (22) 
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Note that the extremization argument h is shifted a.s h + 2{H^r — H"^ Hm) / a"^ h. Although this distribution 
of the vector ^Rt{b — m) is not isotropic, but rather is biased by the matrix \/Rt, the multiphcation by the 
rectangular random matrix H ensures the following approximation under the constraints x — \VRt{b ~ m)\'^/K and 
K = Ke{h^{b — m))/K, where we introduce the auxiliary variables x and k, as follows: 



Tr P{b) exp 

b 



^/R,.S^/Rt{b - m)\^ /a^ + Rc{h^ {b - m)) 
■j dx J d^exp [A' {GstH.H {-x/<J^) + «}] Tr S (\y/Rt{b ^ m)\^ - Kx) 5 {Rc{h\b - m)) - Xk)} ,(23) 



where GgtHrsC^;) is given by Eq. 
Gibbs free energy: 

$(m) 
where ^t{'nn\x) 



Extr 
Extr 

h 



Saddle point evaluation yields the following approximate expression of the 



In Tr ( exp 

b 



KGsiR^s [-^)-Kxx + Mm;x)\+ In Z, 



Rt{b - m)|2 + Re{h\b - m)) 



(24) 



From the Gibbs free energy we obtain a set of equations for estimating m, and using these equations, we construct 
the demodulation algorithm, or the method for finding the minimization argument m for the Gibbs free energy. In 
the following, we summarize the procedure for the minimization of $(m). 

• (Step 0) Initialize variables as x^^^ — 1, m^^'> = 0, h'-^'> — for step t — 0, and set the number of steps as t = 1. 

• (Step 1) For the t-th step, update x a-nd h as 



1 



G' 



• (Step 2) Update m as 



m 



it) = 



where {■)t denotes the expectation 



{f{b))t 



Trb P{b) exp [-x^«)|/R;(b - m)p + Re{h)b)\ f{b) 
Tr5P(b)cxp [-x(*)|Vfi;(6-m)|2 + Re(ft,t6)] 



• (Step 3) Update the number of recursion steps t t + 1. Return to Step 1 unless these variables converge, 
otherwise stop. 

After termination of the above procedure, the transmit signal is estimated as b = argmin^ {1^ ~ m(*)|}. 

The computational cost of Step 1 is 0{KL) and is sufficiently small. In general, the cost of Step 2 is not so small. 
However, we can reduce the cost of Step 2 for the special forms of matrix Rt- For instance, when i?t is a matrix of 
tridiagonal form, as considered in the following, Step 2 can be executed using the transfer matrix method. In such a 
case, the cost is 0{K), which is smaller than the cost of Step 1. 



IV. NUMERICAL EXPERIMENT 



As a simple but nontrivial example, we performed numerical experiments for using Kronecker channel model 
H = \/RiS\/Rt, the correlation matrices i2r and Rt of which are identity and tridiagonal matrices, respectively. 
More precisely, we consider the case in which matrix Rt has nonzero elements only for adjacent antennas, modeled 
hy Rt = I + pR, where the matrix R^k' — Rk'k = 5k(k'-i){k < k') and p is a parameter that represents the strength 
of correlations. For simplicity, we analyze the real channel, and, accordingly, all variables are set to be real. Here, H 
represents the L x X-dimensional i.i.d. Gaussian random matrix A/'(0, 1/L)^^, b = {bk) & is a BPSK-modulated 
transmit signal {bk G {±1}), and rj £ is a normalized real Gaussian-distributed random vector A/'(0, 1). The 
formulation so far for the complex channel can be reconstructed without difficultly for the real channel just by the 
replacement of unitary matrix U with orthogonal matrix O. 
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A. Random-field Ising chain 



Before the analysis of the entire Kroncckcr model, let us evaluate three pieces of mutual information, namely, /i(x)j 
Ii+pO^ Ro{x) Enid Ii+pr{x) Oi that appear in the expression of the mutual information of the entire system Xh, 
as described in Section |IT1 For the BPSK modulation, the mutual information of the one-dimensional channel Ii{x) 
is given by 



hix) 



X~ - I Dz\n (cosh(x + ^/xz)) , 



(25) 



where Dz = exp(— z /2)/v27r. For mutual information //+po^fio(x)i substitution of the tridiagonal form of -R yields 



h+pOTRoix) = Extr |G/+potro (A)+/i(Ax)} 



(26) 



Here, Gj_f_poTiioW = ^(1/2) In (l — (A — l)^/4p^), which is evaluated using the Stieltjes inversion formula for the 
function G{x) (see Reference ^) and the relation pjtt(A) = hme^o(7''-^)~^9Almlndet(ilt — (A — As described 

earlier, the discrepancy between Ii+pr{x) and Ii+po"^ Roix) appears starting from the fourth-order term, which is 
expressed as the tridiagonal form of R, as follows: 



Ii+pr{x) = h+pOTRo{x) 



{-2I'({x)-{2I[{x)fY + ^ 



where 



a(x) = -2 / 7?z{l-tanh2(x + Vx^)}{l-3tanh2(x + ^/^z)} 



(27) 



(28) 



For the tridiagonal form of R, we can evaluate //+pji(x) exactly using the transfer matrix method for the Ising 
chain in random fields. In order to demonstrate how this is accomplished, we transform the mutual information as 
follows: 



I+pR 



(x) = - 



dr' I Tr 
b 2 



K \2ti) 



K/2 



exp 



X In <^ Tr ^ — 
6 2^ V27r 



X 
2 

exp 



r'- ^/l + pRb 



r' - y/TTpR b 
/ DqTilnl Tr TT <j){bk,bk+i\bk,bk+i,r]k) \ 

1 r f 

^ / Tr In <^ Tr TT 



2 



















K-l 



(29) 

where a trivial constant is neglected in the second and the third lines. Note that gauge transformation as 5 ^ r and 
b — > T, where = bkbk and = bk+ibk, and redefinition of r) arc used for simplification of the expression. Matrix 
elements (j){bk,bk+i\bk,bk+i:r]k) and (/)(rfc, r^+i |ffe, 77^) are defined as follows: 



0(6fe,&fe+i|6fc,6fe+i,77fe) = 2*^^? 
(f){Tk,Tk+i\fk,r]k) = -^exp 



\loibk - bk) 4- li{bk+i - bk+i)\ + y/xVk (loh + hh+i) 

\lo(Tk - 1) + hTk{Tk+l - 1)1^ + ^/XVk {loTk + hfkTk+l) 



(30) 



where Iq and li are real constants that are obtained by Cholesky decomposition, i.e., / -t- pR = AA^, where Kkk = 
^0, A(fe+i)fc = /i, and zero otherwise. These constants satisfy = 1, ^o^i = and Iq > li. The matrix element of 

Eq. (|30p corresponds to the Boltzmann weight of the Ising chain coupled with bimodal (r) and Gaussian (ry) random 
fields, and consequently the bit error rate (BER) for the demodulation bk — argmaxf,^ {Trb\{,^ P(b|r)} is evaluated 
analytically for the random field Ising chain. 

Several methods of analysis have been developed for handling the random field Ising chain [l^, [23|, and in the 
present study, we use the technique of belief propagation, which is equivalent to the transfer matrix method. After 
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parameterization of the belief from site k to site fc + 1 by cavity field h^k+i as /i(Tfc+i) = 6^'^''+'^'^''+'^ /2 cosh{h^k+i), 
the propagation process from to h^^+i is written as 

/ . X 1 .1 fT.r,r, + ^'l^iTk,Tk+l\fkr]k)exp{h^kTk)Tk+l\ 

h^k+1 = arctanh } fi{Tk+i)Tk+i = arctanh — = — p — r — 

J \ Y.r,T^ + ^nTk,T-k+l\TkVk)exp{h^kTk) J 

= Xlf + XPn + y/xhTkVk ~ ffcarctanh(tanh(/9x) tanh(/i^fe + x^o + XP^fc + ^loVk))- (31) 
Similarly, for the opposite direction, denoted by hk^, 

h ~ ^ ^fT.r,Tk+,(f>{Tk,Tk+l\fkr]k)exp{hk+l^Tk+l)Tk\ 

''^ \ T,T,T,+,HTk,Tk+i\fkrik)exp{hk+i^Tk+i) J 

= xIq + XRTk + y/xloVk - ffearctanh(tanh(px) ta.-nh{hk+i^ + x^i + XPT'k + y/xhfkVk))- (32) 

The stationary distributions of beliefs for the numerically increasing and decreasing directions, denoted by 7r+ and 
7r_, respectively, satisfy the following conditions: 

TT+{h^k+i) = j 'K+{h^k)dh^k J Dt] ^d^h^k+i-ixli+XPT + y/xhfr]) 

+f arctanh (tanh(yox) tanh^h^k + X^o + XP^ + Vx^ov)) ^ > (33) 
■n-{hk^) = j ■n-{hk+i^)dhk+i^ / X! ^/ife^ - (x^o + XP^" + Vx^o??) 

+f arctanh (tanh(yox) tanh(/n:+i^ + x^i + XP^" + VxhTf])) ^ , (34) 

where / 7r+(/i_»)(i/i^ = J Tr^(h^)dh^ = 1. The functions Tr^{h^k+i) and 7r_(/ifc^) can be obtained numeri- 
cally by the Monte Carlo method for a one-dimensional system. The bit error rate Pb is represented by Pb = 
J Tr+{h^)dh^ J Tr_{h^)dh^ {1 - sgn{h^ + h^)) /2. 

In order to investigate the performance of communication channels, conditional entropy h{b\r) — h{b) —Iff, where 
h{b) denotes entropy, is a favorable measure, because h{b\r) decreases to zero under smaller noise power. Figures 
[1] and show the exact conditional entropy hi^pR{b\r) estimated by the Monte Carlo method, the approximate 
entropy ^f+po^RoC'l''') obtained by matrix integration, and the entropy of the i.i.d. channel hi[b\r). In both 
graphs, the entropy obtained by matrix integration, i.e., hi^pQT Ro{b\r), does not exceed the entropy obtained by 
exact evaluation hi+pii{b\r)^ which implies the inequality Ii+pR(b,r) < Ii+po'^Roib,'''), which is given up to the 
fourth order of perturbation in Section [TTl The analysis indicates that the deviation of the approximate result from 
the exact result depends on the signal-to-noise ratio (SNR) and the correlation parameter p. For a small SNR and 
small correlation (= small p) the deviation is small, which means that the entropy, ^f+po^fto(^l'") obtained by matrix 
integration gives a good approximation of the exact entropy, /if+pfi(6|r), while the deviation becomes greater in the 
case of a large SNR or large correlation. 

As we have discussed earlier, the approximate evaluation with matrix integration is useful because this simplifies 
the analysis. However, as shown by the numerical results for conditional entropy, this method is only valid when the 
correlation is small. 



B. Kronecker channel 



Next, we examine whether the proposed demodulation algorithm is practical. In Figs. [3] and [H the results of 
demodulation for the Kronecker channels are depicted. These figures show two BER curves, namely, the BER obtained 
by exact analysis of the model Rt = I + pR with tridiagonal R, as described in the previous subsection, and the BER 
obtained by the correlation matrix multiplied by an arbitrary orthogonal matrix, Rt — I + pO^ RO and averaged 
over the orthogonal matrix O. For the latter model, matrix integration analysis can be applied due to multiplication 
of the orthogonal matrix. We have proposed an appropriate demodulation algorithm for each evaluation. For the 
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FIG. 2: Conditional entropy hi+pii{b\r) for chain-like system (2). BER vs. correlation parameter p ((a) SNR — 2 dB, (b) 
SNR — 6 dB). The solid lines show the results obtained by exact analysis, the broken lines show the results obtained by matrix 
integration with correction from the fourth order, and the dotted lines show the results obtained by matrix integration without 
correction. 



former model, without the orthogonal matrix multiplication, we can use the belief propagation algorithm proposed in 
the previous subsection. For the latter model, with orthogonal matrix multiplication, the demodulation algorithm we 
proposed in based on matrix integration and the Thouless- Anderson-Palmer method [2l|, [g^, HH is applicable. 

The results of demodulation for each case show good agreement with the results obtained by replica analysis. For 
large noise power and large p (Figs. [D^b) and[4l^b)), the BER of the model without orthogonal-matrix multiplication 
becomes larger than that with orthogonal-matrix multiplication, which reflects the discrepancy of mutual information 
from higher-order values of p, as mentioned in Section [TTl Therefore, in designing the demodulation algorithm for 
such correlated channels, appropriate treatment of the correlation matrix should be taken into consideration. We also 
examined the convergence speed of the algorithm. The proposed algorithm in the previous subsection requires the 
0{K^) matrix operation and the 0{K) belief propagation process in each step, and convergence of this algorithm 
requires dozens of iterations. Therefore, we conclude that this algorithm is computationally feasible. 



V. SUMMARY 



In the present paper, we proposed a performance analysis scheme and a demodulation algorithm for the Kronecker 
channel model in MIMO wireless communication systems. For a more exact evaluation than that of our previous 
study using the matrix integration formula, we demonstrated that two separated manipulation steps for the product 
form of the channel matrix, i.e., Gaussian integration for the channel matrix and an appropriate scheme for dealing 
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FIG. 3: Performance by the replica analysis and the result of demodulation experiment (1). BER vs. SNR ((a) p — 0.2, (b) 
p — 0.5). We set the parameters K — 4, 400 and L = 4, 000 and take the average of the results over 128 samples with various 
input signals b, matrices H, and noises rj. We also varied the orthogonal matrix O for the model with orthogonal-matrix 
multiplication. The lines depict the results of the replica analysis for two models, the chain-like model Rt — I + pR with 
tridiagonal R (solid) and the model with orthogonal-matrix multiplication Jtt = / + pO^ RO (dotted). The symbols in the 
figure denote the results of demodulations, namely, demodulation for the chain-like model (*), demodulation for the model with 
orthogonal- matrix multiplication (x), inappropriate choice of the demodulation algorithm, i.e., the demodulation algorithm for 
the orthogonal- matrix multiplication model applied to the chain- like model (+). 
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P P 

FIG. 4: Performance by the replica analysis and the results of demodulation experiment (2). BER vs. correlation parameter 
p ((a) SNR = 5 dB, (b) SNR = 7 dB). We set the parameter K = 4,400 and L = 4,000 and take the average over 1,024 
samples for various values of 6,H,t7, and O. The lines depict the results of the replica analysis for two models, the chain-like 
model Rt = I + pR with tridiagonal R (solid) and the model with orthogonal-matrix multiplication Rt = I + pO^RO 
(dotted). The symbols in the figure denote the results of demodulations, namely, demodulation for the chain-like model (*), 
demodulation for the model with orthogonal-matrix multiplication (x), inappropriate choice of the demodulation algorithm, 
i.e., the demodulation algorithm for the orthogonal-matrix multiplication model applied to the chain-like model (-I-). 

with transmitter correlation, are important for the correlated MIMO system. The numerical result for the tridiagonal 
correlation matrix model shows that the proposed scheme and algorithm are useful for performance analysis and for 
the construction of a practical demodulation algorithm. 

APPENDIX: PERTURBATIVE EXPANSION OF FREE ENERGY 

In the appendix, we derive the perturbative expressions of free energies from two evaluations, namely, matrix 
integration and exact analysis. Remember that the mutual information for a MIMO system can be obtained by 
evaluating Eq. ([6]). For convenience, we introduce the constant x iiito the channel definition (see also Eqs. (fT5|) and 
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((T6)) ). as follows: 

r = ^Hb + ri, (A.l) 

where we set a = 1 in Eq. ([T]). Taking QPSK modulation into account, let us assume the probability distribution 
P{bk), which satisfies the following conditions: 

• 1. P{bk) can be factorized into the same distributions for the real and the imaginary parts: P{bk) = 
P(Re(bk))P(Im(6fe)). 

• 2. Reflection symmetry: Fix) = P{—x). 

• 3. Signal power condition: Y.b, P{^c{hk))^c{hkf = J2b, ^^(Ini(fefe))Im(6fc)' = 1/2. 

From these conditions, we have J^b, Pih)bf-' = 0, J^b, P(bk)\bk\''^bk ^ ,Eb, P{h)\bk\^ = 1 {I E N), and so on. 
Quadrature-phase shift keying modulation is included in this case. As mentioned in the text, we also assume that 
the cross-correlation matrix H'^H can be written as H^H = I + pR with zero diagonal elements of R, R^k — and 
convergence of the eigenvalue distribution of the cross-correlation matrix for K ^ oo. 



1. Expansion of free energy via matrix integration 

As shown in Section [TTl free energy is obtained via matrix integration as follows: 

Ii+pU^Ruix) = Extr iliilx) + il- Gi+pr{0] ■ (A.2) 

The function Gi+pii{z) can be decomposed to obtain 

Gi+pr{z) = z + Gr{pz). (A.3) 
Let us define A" = Tr(ii")/7^. The function Gr{z) can be expressed in terms of A", as follows: 

Gr{z) = I'^z^ + I'^z^ + 7 - 2A^') + • • • , (A.4) 



2 3 4 

from the formula G(z) = dx{f{x) - x'^) s.t. x = J d\p{X){f{x) - \)-\ Note that A = for and G(0) = 
Substitution into mutual information after redefinition of ^ yields the following: 

Ii+pUiRuix) = Extrf//(x + ex) + Ce-Gfi(pe)| 

= Extr I + hl'iix) + ^^Iiix) + • • • 



(A.5) 



From the saddle point conditions with respect to we have 

I - pG'r {pO = pG'r i-pxi'iix) + cp^) = pG'r {-pxI'Ax)) + cp^G^ (-px/;(x)) + O(p^) 

= -p^yxl'iix) + P'A^X'{/;(X)}' - P\~^ - 2X^\^{I'j{x)Y + cp^A^ + 0(p'^), (A.6) 

for p up to the fourth order (c is an 0(1) constant). Substituting these equations into the original expression for the 
mutual information, we obtain 

h+pu^Huix) = //(x)-^A^{/;(x)}^ + ^A^{/;(x)r-^(A^-2A^^){/;(x)r 

+^A^'{/z (X)}{/;(X)}' + 0{p'). (A.7) 
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2. Expansion of rigorous free energy 



Before studying perturbative expansion of the complex channel, let us start with the real channel for which mutual 
communication is given by 



1 



X 



r-\ Rtb 



InTr 

b 



K/2 



X 



exp 



f - \ Rtb 



;in(^)-; 



(A., 



where all variables and matrices are real and denoted with tilde for discrimination between the real and the complex 
channels in this subsection. By substituting Rt — I + pR into the above equation, where i? is a symmetric matrix 
with zero diagonal elements, expanding with respect to p, and then performing some algebraic manipulation, we 
obtain the following: 



I+pR^XJ 



k=l 



{PX? 1 
3! 2K 

{pxf 1 
4! 2K 



RiijiRi2j2^i3j3 ([{biibi2bi3)c{bjibj^bj^)c] - A{biibj2)c{bi J) j^)c{h J) ji)c\j 



E 



Riiji Ri2j2 Risjs Riiji 



ni2l3«4jlj2j3j4 

X (^{KKKK)c{bhbjJ>jJ)j^)c] - 12[(6,i6j,6j3)c(&ji&j,6jjc(fe»3^j4>c] + 6[(&,i&j,)c(fe»2^j3>c(^»3^i4>c(&«4^Ji>c] 



= ir\x) 

[pxY 



K 



2\ 13 



{pxf E.,-4 [^^4^^]2 



K 



En4 



if X 

A3 [(6^ 



6 K ' ' 48 X 



K 



if 



2\2ir/l,2\ 12 



^^^^ V[(62)c]^ + ^A^[(&^)c]^ - ^A^[(6'> 
6 8 



K 



K 



2\ 12\2 



+ 0{p^). (A.9) 



Here, = Tr^ P(6)/(b) exp(-x|r - 6| 2/2)/ Tr^ P(b) exp(-x|r - 6| 2/2), and (/(b))c is its cumulant. In addition, 

{fib)) is a function of f , and is always accompanied by the average [-F'('r)] = (x/27r)^/2 J df. P{b)F{f) exp(— xl'r — 
6p/2). where & without a subscript denoted the signal at an arbitrary antenna. The expression in the third line of 
Eq. (|A.9[) is obtained by repeatedly performing integration by parts with respect to f. We can also show that 
= 2/^'-°'^'(x) and [{lYA = -2Ii""'\x), from which we have 

24Tpr(x) = 2/r'(x) - ^A^{2//^^i(x)}' + ^A^{2/r^'(x)}' - ^F{2//-'(x)}* 



(PX)M2E.(^')? 



(-2//--i(x) - {2i'r\x)r){2i'r\x)r 

[{b%?y 



K 



ij ^ij 

K 



(-2/;'--i(x) - {2/r^'(x)}')' + 



'roal/,,M 2\2 



6 



Oip' 



(A.IO) 
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Next, we convert this result into the result for a complex channel, the mutual information of which is given by Eq. 
P^ . We can easily show that the complex channel described by r = ^yxVRtb + 2; is equivalent to the real channel 

of double size f = ^/xy/Rtb + i, where 



For such a system, we can show that //(x) — 2/}° (x) and Ir^{x) = 2/^ (x)- Since the eigenvalue distributions of 

Rt = I + pR and the corresponding are the same, from the relationship between the real and complex channels, 
we have 



ii+pRix) = 2j-;^(x) 



+0(p5). (A.12) 
The function C(x) is given by 

dx) ^ 2{xM^ J drTrP(b)(Re(6)4 + Im(&)4)^-Pexp(-x|r - bp), 
where (/(b))™P = Tr P(b)/(b) exp(-x|r - bp)/ Tr P(6) exp(-x|r - bp), (A.13) 

b b 

and the subscript of the angular bracket c denotes the cumulant. Substituting the definitions of and A^, the 
discrepancy between the two results is obtained as 

ii+pR{x)~ii+,WRuix)^ E {(^')»» - ^^}' {-iiM I'lix)'} {i'Ax)r 



4 ii: 



^{Re(i?.,)4 + Im(i?,,)n {-/;'(;^)-/;(;,)2}^ + ^ +0(p5)^ (A.14) 



that is, the dominant term of the discrepancy is of the order p . The factor —Ij{x) ~ {I'lix)} = ^'^Ij'^ ix) ^ 
|2j^rcai^^^|2 _ [(fe^)^] — [{b^)c]^ is uonncgative, and the inequality Ii+pu^ Ruix) ^ Ii+pr{x) holds up to the fourth 
order. 
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